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FAST  ALGORITHMS  FOR  PARTIAL  DIFFERENTIAL 
EQUATIONS  ON  ADVANCED  COMPUTERS 


Final  Report  AFOSR  86-0061 
Thomas  A.  Manteuffel 

Introduction. 

This  report  will  cover  four  areas  of  research.  The  first  area  is  an  analysis  of 
preconditionings.  We  have  examined  the  properties  that  a  preconditioning  must 
possess  in  order  to  yield  an  iterative  method  with  convergence  properties  indepen¬ 
dent  of  the  discretization  parameters.  The  second  area  is  an  attempt  to  provide 
a  framework  for  the  construction  of  conjugate  gradient  methods.  The  third  area 
of  research  is  the  analysis  of  the  equations  governing  the  transport  of  neutrally 
charged  particles  and  the  construction  of  iterative  methods  for  the  solution  of  dis¬ 
crete  transport  equations.  The  final  area  of  research  is  supraconvergence  which 
deals  with  accurately  approximating  the  solution  of  partial  differential  equations  on 
highly  irregular  meshes. 

The  first  topic  is  very  theoretical  in  nature  but  will  have  important  implications 
to  the  practical  choice  of  preconditionings.  The  second  area  is  also  theoretical  in 
nature,  but  leads  immediately  to  the  development  of  new  iterative  methods.  The 
third  area  builds  on  the  analysis  of  the  continuous  operators  to  yield  insight  into 
the  behavior  of  discrete  approximations.  The  final  goal  of  this  third  project  is  the 
development  of  efficient  numerical  methods  for  the  solution  of  discrete  transport 
equations  on  massively  parallel  machines.  The  final  area  of  research  will  lead  to 
more  efficient  solution  of  differential  equations  on  irregular  meshes.  v  p  j 

l.  Analysis  of  Preconditionings 

Uniformly  elliptic  partial  differential  equations  in  R2 ,  denoted  here  by 


yield  discrete  linear  systems  of 


--  /, 


f. 


with  condition  number 

■■  '  0(h~2), 

where  h  is  a  mesh  parameter.  I'*r.i’:\e  methods  for  the  solution  of  the  discrete 
equations  yield  error  bounds  after  m  of  the  type 


u.;j 


]e  M 
1 1 


89 


•? 


where  k  =  Cond{Ah)  or  Cond[Ah) »  depending  on  the  implementation. 

Preconditioning  can  be  viewed  as  a  linear  rearrangement  of  the  equations  in 
an  attempt  to  reduce  the  condition  number.  The  discrete  equation  is  replaced  by 

ChAhu  =  Chl 

where  Ch  is  a  nonsingular  linear  process.  Classical  splitting  methods  may  be  viewed 
from  this  perspective.  For  example,  S.O.R.  with  optimal  parameter  (applied  to 
appropriate  systems)  yields  a  preconditioning  for  which 

Cond(ChAh)  =  0(h~x). 

Likewise,  incomplete  factorization  techniques  yield  preconditionings  with  this  same 
property.  Although  of  the  same  asymptotic  order  as  S.O.R.,  incomplete  factoriza¬ 
tions  perform  better  on  a  wide  class  of  problems  because  the  absolute  magnitude 
of  the  condition  number  is  smaller. 

It  was  shown  by  D’Yakanov  DYl’,  ;DY2j  and  Gunn  [GUI]  that  for  a  special 
class  of  self-adjoint  uniformly  elliptic  operators  it  is  possible  to  construct  Bh  such 
that 

Cond(B~l  Ah)  =  0(1), 

that  is,  the  condition  number  of  the  preconditioned  system  is  independent  of  the 
mesh  parameter.  Here  Bh  is  the  discretization  of  a  ‘nearby’  uniformly  elliptic  op¬ 
erator,  B.  The  error  in  the  original  problem  can  be  reduced  by  a  given  factor  in  a 
fixed  multiple  of  the  work  required  to  solve  systems  of  the  type 

Bhu  =  f. 


The  multiple  is  independent  of  the  mesh  but  depends  upon  how  ‘close’  A  is  to  B. 

This  result  was  exploited  by  Concus  and  Golub  [C01]  and  Widlund  [WIl]  who 
noticed  that  if  A  and  B  live  on  a  rectangle  and  if  B  is  separable,  then  the  equation 
involving  Bh  can  be  solved  using  fast  direct  methods  (cf.  [SWl]).  Practice  has 
shown  that  these  methods  are  of  little  use  unless  A  is  very  ‘close’  to  a  separable  B. 
The  difficulty  is  that  while  Condi  B~  1 Ah )  is  independent  of  the  mesh,  for  reasonable 
values  of  the  mesh  parameter  h,  the  magnitude  of  the  condition  is  larger  than  that 
for  other  splittings,  for  example,  incomplete  factorization. 

Multigrid  may  ?dso  be  viewed  as  a  preconditioning  where  Ch  represents  the 
process  of  cycling  through  the  grids.  For  a  wide  class  of  problems,  Cond(ChAh )  is 
not  only  independent  of  h,  but  very  small. 

Results.  The  methods  of  D’Yakanov  and  Gunn  and  multigrid  methods  satisfy  the 
same  asymptotic  bounds  on  the  condition  number.  Yet,  in  practice,  their  perfor¬ 
mance  can  differ  greatly.  An  attempt  to  develop  a  theory  that  encompasses  these 
methods  and  yields  more  precise  measurements  of  the  practicality  of  precondition¬ 
ing  was  presented  in  a  report  by  Faber,  Manteuffel,  and  Farter  FA4].  This  work 
developed  the  concept  of  equivalent  operators  on  Hilbert  spaces.  Two  operators  A 

Ulit 

I  \h 


a 

*  □ 

vd  □ 


■t 

..■b’tty  Codes 

.ail  and /or 
Special 


and  B  are  said  to  be  equivalent  in  norm  on  the  subset  D  of  the  Hilbert  space  H  if 
the  ratio 

!!ax|!/||Bx|| 

is  bounded  above  and  below  for  every  element  xeD.  Equivalence  in  norm  yields 
Cond(AB~  1 )  bounded  on  an  appropriate  domain.  For  preconditioned  systems  it 
is  more  appropriate  to  examine  Cond{B~  1  A)  or  to  show  that  A~  1  and  B~  1  are 
equivalent  in  norm.  An  important  result  in  this  report  is  that  if  Ah  and  Bh  are  a 
sequence  of  discrete  approximations  to  A  and  B,  then  Cond(B^  1  Ah)  is  bounded 
independent  of  h  only  if  Cond(B~  1  A)  is  bounded.  Thus,  the  construction  of  pre- 
conditionings  independent  of  the  mesh  must  be  based  upon  equivalent  operators. 

This  report  was  the  basis  for  two  papers  [MA5],  [FA5]  (attached).  In  the  first, 
uniformly  elliptic  operators  were  examined.  It  was  shown  that  if  fi  is  a  sufficiently 
smooth  region  in  R2 ,  then  any  two  uniformly  elliptic  operators  with  the  same  ho¬ 
mogeneous  boundary  conditions  are  equivalent  on  a  dense  subset  of  L2  ( f2).  One 
interesting  corollary  of  this  result  is  that  if  Ah  and  Bh  are  finite  element  discretiza¬ 
tions  of  uniformly  elliptic  operators  using  elements  smooth  enough  to  be  in  the 
domain  of  A  and  B,  then  A h  and  Bh  are  uniformly  equivalent;  that  is,  they  are 
equivalent  and  the  bounds  are  independent  of  h.  This  implies  that  Cond[AhB^ 
will  be  independent  of  the  mesh.  This  will  lead  to  bounds  on  the  residual  rm  =  Aem 
of  the  type 

r"  im 


where  k  =  Cond(AhBh  1). 

To  obtain  similar  bounds  on  the  error  em  it  is  necessary  to  examine  Cond 
( B~lA ).  It  was  discovered  that  Cond{B~  1  A)  is  bounded  if  and  only  if  A*  and  B* , 
the  L2-adjoints  of  A  and  B,  have  the  same  homogeneous  boundary  conditions.  It  is 
easy  to  construct  examples  of  operators  A  and  B  where  Cond(AB~  l)  is  bounded 
while  Cond(B~  1  A)  is  not  and  conversely.  This  discovery  negates  the  conventional 
wisdom  of  choosing  B  to  be  the  leading  part  of  A  with  the  same  boundary  condi¬ 
tions. 

The  second  paper,  [FAS’,  extends  the  theory  in  [FA4]  to  cover  the  situation 
where  A  and  B  are  operators  from,  a  Hilbert  space  W  to  another  Hilbert  space  V. 
This  paper  establishes  in  a  more  general  context  the  same  results  as  [MA5l.  Finally, 
this  paper  examines  the  uniform  <-q  valence  of  finite  element  approximations.  Sim¬ 
ilar  results  are  established  for  so::..-  -Mndard  finite  difference  approximations  when 
fi  is  a  rectangle. 

These  two  papers  also  exa::  •  >■  ■  q  oivalence  of  weak  forms  of  the  operators 
A  and  B.  Here,  however,  the  //.  ::::  laces  the  L2-norm  in  the  definition  of 

equivalence.  It  was  discovered  i. ii.it  »  :  B‘l A)  (the  Bj  condition  number)  and 

CWl(AB_1)  are  both  bounded  if  and  only  if  .4  and  B  have  Dirichlet  boundary 
conditions  on  exactly  the  same  portions  of  the  boundary  of  fl.  Similar  results  hold 
for  their  discrete  counterparts. 
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Future  Research.  These  results  imply  that  under  very  general  hypotheses  precon¬ 
ditioning  the  discrete  approximation  of  a  uniformly  elliptic  operator  by  the  discrete 
approximation  of  any  other  uniformly  elliptic  operator  with  the  appropriate  bound¬ 
ary  conditions  will  yield  a  condition  number  independent  of  the  mesh.  This  result 
demonstrates  that  a  mesh  independent  bound  alone  is  not  sufficient  to  yield  a  good 
preconditioning.  A  measure  of  the  distance  between  operators  has  been  suggested 
[FA5j  but  is  of  little  practical  value.  Conversely,  these  results  imply  that  pre¬ 
conditioning  strategies  that  yield  mesh  independent  condition  numbers  most  likely 
approximate  the  inverse  of  some  operator  equivalent  to  A. 

A  long  list  of  interesting  questions  arises,  some  of  them  theoretical  and  some 
with  more  practical  application.  The  above  results  depend  upon  H2  regularity  (or 
H x  regularity  for  the  weak  forms)  and  uniform  ellipticity.  One  theoretical  question 
is  whether  these  results  hold  in  the  absence  of  regularity.  Another  question  is  how 
to  choose  the  operator  B  so  that  B~  1 A  has  its  spectrum  in  the  right-half  plane,  a 
potentially  important  issue  for  nonself-adjoint  operators.  This  would  determine  the 
type  of  iteration  implemented  and  whether 

k  =  Cond(B ^  1  Ah )  or  Cond(B ^  1  Ah )  ’ 
in  the  error  bound  above. 

Of  more  practical  interest  would  be  the  development  of  a  measure  of  the  dis¬ 
tance  between  two  elliptic  operators  that  could  be  computed  or  at  least  well  ap¬ 
proximated  a  priori.  This  would  yield  a  good  approximation  to  Cond(B~  1  Ah)  and 
allow  one  to  accurately  predict  the  work  required  to  solve  the  discrete  system. 

Even  more  fruitful  would  be  the  examination  of  multigrid  methods  in  this  con¬ 
text.  Recent  work  by  Parter  [PAl]  on  singularly  perturbed  Helmholtz  type  equa¬ 
tions  has  shown  that  in  some  cases  one  can  achieve  a  mesh  independent  condition 
number  using  partial  multigrid  methods.  These  methods  keep  the  coarsest  grid 
very  fine  and  never  completely  solve  the  coarsest  mesh.  This  leads  one  to  believe 
that  similar  methods  may  be  applied  to  other  problems.  The  marriage  of  multigrid 
with  iterative  techniques  may  guarantee  a  mesh  independent  condition  number  for 
a  much  wider  class  of  problems  than  those  for  which  multigrid  has  been  applied  to 
date. 


2.  Con'jd;  gradient  Methods. 

Since  the  introduction  of  the  mr’.jugate  gradient  method  in  1952  by  Hentenes 
and  Stiefel  [HEl]  it  has  become  a  pop  Far  method  for  the  solution  of  symmetric  pos¬ 
itive  definite  linear  systems.  During  ‘hr  36  years  since  its  introduction,  many  vari¬ 
ants  have  been  developed  includinu  preconditioned  conjugate  gradient  method 
and  the  conjugate  residual  metho<i 

Efficient  generalizations  of  the  conjugate  gradient  method  to  nonself-adjoint 
linear  systems,  other  than  forming  the  normal  equations,  were  sought  but  never 
found.  In  the  early  80’s  work  by  Faber  and  Manteuffel  laid  to  rest  this  question 
by  establishing  both  necessary  and  sufficient  conditions  for  the  existence  of  a  short 
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term  conjugate  gradient  method  [FAl],  [FA2] .  In  brief,  given  a  linear  system  A 
and  a  preconditioning  C  it  is  possible  to  construct  a  three  term  conjugate  gradient 
method  optimal  in  the  innerproduct  norm  ||x|]|  =  (Bx,x)  if  and  only  if  i)  CA 
has  fewer  than  three  distinct  eigenvalues,  or  ii)  CA  is  5-normal  (1).  The  second 
criteria  is  an  extension  of  the  standard  definition  of  normality.  Here  normality  is 
measured  with  respect  to  the  inner  product  { B- ,  •).  The  criteria  is  equivalent  to  the 
expression 

BA  =  aB  -f  (3 A'  B 

where  A*  is  the  standard  adjoint.  This  work  showed  that  the  conjugate  gradient 
method  could  be  extended,  but  only  to  a  relatively  small  class  of  nonself-adjoint 
linear  systems. 

Results.  The  paper  by  Ashby,  Manteuffel  and  Saylor  [ AS  1  ]  (attached)  is  an  at¬ 
tempt  to  use  the  theory  developed  in  fFAl]  and  [FA2]  to  provide  a  structure  that 
includes  all  conjugate  gradient  methods.  This  structure  not  only  shows  how  many 
variants  of  the  conjugate  gradient  method  are  really  different  implementations  of 
the  same  method,  but  also  unambiguously  determines  the  set  of  linear  systems  for 
which  the  method  is  applicable.  Finally,  this  theory  leads  to  the  easy  development 
of  new  methods.  In  fact,  this  paper  introduces  several  new  variants  of  the  conjugate 
gradient  method. 

Future  Research.  The  same  principles  have  been  applied  to  a  much  broader  class 
of  methods  known  as  Projection  Methods  which  include  the  conjugate  gradient 
methods  and  are  applicable  to  nonself-adjoint  linear  systems.  Recent  work  by  Jou- 
bert  and  Manteuffel  [JOl]  has  established  the  beginnings  of  a  theoretical  framework 
for  these  methods.  Future  research  will  include  expanding  these  results  and  exam¬ 
ining  several  new  methods  exposed  by  this  theory. 

3.  Multigrid  Methods  for  Transport  Theory 

Mathematical  models  of  radiation  transport  in  optically  dense  materials  are  a 
fundamental  part  of  many  weapons  codes.  Similar  problems  arise  in  nuclear  reactor 
design  and  analysis,  satellite  electronics  shielding,  radiation  detector  design  and 
dosimetry,  radiation  therapy  planning,  and  radiation  effects  on  materials.  Although 
two  and  three  dimensional  anisotropic  models  are  of  great  interest,  even  the  one 
dimensional  model  of  transport  in  slab  geometry  presents  formidable  difficulties. 
This  is  because  the  equations  governing  transport  become  nearly  singular  on  a 
subspace  of  large  dimension  in  optically  dense  material,  that  is,  when  the  mean 
free  path  between  particle  collisions  is  small  compared  to  the  physical  dimensions. 
This  is  further  complicated  in  that  the  subspace  contains  the  functions  of  lowest 
frequency  in  space. 

To  be  more  precise,  consider  the  equation  governing  neutron  transport  with 
isotropic  scattering  in  a  slab  geometry, 

d  7  f 1 

+  u(m,x)  =  -  /  u(fj.}x)dn  + q. 
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Here,  u  is  the  particle  density  at  spatial  position  x  moving  in  direction  characterized 
by  angle  #,  where  p  =  cos#,  and  q  is  the  source  term.  This  equation  is  to  hold  for 
x e(  —  a,a)  with  boundary  conditions 


u(p,—  a)  —  6j  (p.)  0  <  p  <  1 

u(p,a)  =  b2  (/n)  -  1  <  p  <  0. 

Here,  7  is  a  positive  constant  with  7  <  1.  The  space-like  parameter  x  is  scaled 
in  terms  of  the  number  of  mean  free  paths.  Problems  of  interest  generally  involve 
7  =  1  and  very  large  a. 

The  problem  can  be  reformulated  by  defining  the  operators 


L-  = 


•dp 


and 


yielding 


/, 


Hu  u  =  -)Lu  -r  q. 


Inverting  the  transport  operator  Htl  yields 


u  =  H~  1 7LU  +  H~lq. 

Applying  L  to  this  equation  and  making  the  substitutions 


«^  =  Zu,  /  =  LH~  1  q, 

yields 

*  -  = 

The  operator  K  =  ( LH~  *7)  can  be  shown  to  be  a  compact  self-adjoint  oper¬ 
ator  on  L2{— a,  a)  with  \\K\\  <  1  n  f.  VVT2j).  Notice  that  applying  L  eliminates 
the  dependence  on  the  angle  p.  Tin*  equation  above  can  be  viewed  as  an  integral 
equation  of  the  second  kind  on  /._  a.  u). 

For  very  large  u,  a  singular  value  decomposition  of  K  yields  singular  vectors 
that  are  asymptotic  to  Fourier  vec.irs. 

I  ">*  '  \  )  k  odd 

Vk  1  '  i  k  even 

where 

k  k 

r)k  =  kn  -  0(  -  ),  for  -  <<  1, 
a  a 
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and  singular  values 


*-!-*<£)* +o((;)*). 

where  1/4  +  —■  <  6k  <  1/4  -f  ~  [FA3]  (attached). 

It  is  often  the  case  that  the  solution  can  be  well  described  by  frequencies  for 
which  -  <<  1.  This  implies  that  the  operator  I  —  K  is  nearly  singular  on  frequencies 
of  interest. 

The  discrete  analogue  of  this  equation  retains  many  of  these  properties  and 
introduces  greater  difficulties.  The  operators  H M  and  L  are  discretized  on  a  ( x,/u ) 
grid.  A  popular  technique,  known  as  the  SN  or  discrete  ordinates  method  [CAl], 
is  to  approximate  Hu  by  a  difference  formla  in  space  and  to  approximate  L  by 
a  quadrature  formula  in  n.  Thus,  the  discrete  Kh  can  be  applied  by  computing 
(Hilx)~l4>h  for  a  discrete  set  of  angles  /z,  and  applying  a  quadrature  rule  to  form 
Khih. 

In  the  past,  these  problems  were  solved  by  the  simple  splitting  iteration  known 
as  the  source  iteration, 

1  =  +  fH. 

For  large  a,  this  iteration  will  converge  very  slowly.  Recently,  an  observation  that 
slow  convergence  occurs  in  problems  that  are  basically  diffusive  has  led  to  an  algo¬ 
rithm  called  Diffusion  Synthetic  Acceleration  (DSA)  (cf.  Larson  [LAl],  [LOlj), 

(■ I~Gh )<brl  ={Kk-Gh)<t>ih+fh, 

where  Gh  —  1  and  Th  is  the  discrete  form  of  the  diffusion  operator 


Tu  —  — 


1  d2u 
3  3X7 


-+-  u 


with  u  satisfying  appropriate  boundary  conditions.  The  DSA  iteration  can  be  re¬ 
arranged  to  yield 


#+1  =  (Th  -  I)~'-{ThKh  -I)<ph  +  (Th  -  I)~lThf . 

Thus,  the  iteration  requires  only  the  solution  of  a  two-point  boundary  value  problem 
and  application  of  the  operator  K  at  each  step. 

Results.  Recent  work  by  Faber  and  Manteuffel  [FA3]  (attached)  includes  an  anal¬ 
ysis  of  the  singular  values  and  singular  vectors  of  the  operator  K  in  the  limit  of 
large  a.  These  results  are  used  to  that  the  DSA  iteration  is  successful  because 
the  first  few  singular  vectors  and  singular  values  of  the  selfadjoint  compact  operator 
K  are  well  approximated  by  the  first  few  -dngular  vectors  and  singular  values  of  the 
compact  operator  G.  Put  in  this  perspective.  (/  —  G)  is  a  good  preconditioning  for 
(I  —  K).  Further,  simple  iterative  methods  can  be  used  to  accelerate  convergence  of 
the  DSA  iteration.  In  fact,  numerical  experiments  show  that  a  conjugate  gradient 
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iteration  with  DSA  as  a  preconditioner  can  reduce  the  overall  work  by  a  factor  of 
three  even  in  the  most  simple  cases. 

The  DSA  algorithm  is  extremely  effective  for  problems  with  isotropic  or  mildly 
anisotropic  scattering,  but  it  is  ineffective  for  problems  with  highly  anisotropic 
scattering,  such  as  charged  particle  problems  [MOlj.  Furthermore,  certain  practical 
difficulties  often  arise  in  applying  the  DSA  method  with  high-order-accurate  spatial 
differencing  schemes.  Specifically,  stability  considerations  require  that  the  diffusion 
operator  used  in  the  DSA  scheme  be  derived  directly  from  the  spatially-differenced 
SN  equations,  and  non-standard  forms  of  the  diffusion  equation  that  are  very  costly 
to  solve  can  be  obtained  [ALl], 

In  principal,  multigrid  methods  are  well  suited  to  the  solution  of  integral  equa¬ 
tions  of  the  second  kind.  Because  K  is  compact,  for  any  e  >  0  all  but  a  finite  number 
of  the  singular  values  of  K  are  less  than  e.  Thus,  I  —  K  looks  like  the  identity  plus 
a  finite  rank  operator.  In  addition,  the  lowest  frequency  singular  vectors  are  those 
associated  with  the  largest  singular  values  of  K .  Potentially,  these  vectors  could  be 
well  approximated  on  a  coarse  grid. 

For  large  o,  however,  all  frequencies  of  interest  may  be  associated  with  singular 
values  of  K  very  Hose  to  unity.  This  difficulty  is  compounded  in  forming  the  discrete 
analogue.  If  the  spatial  mesh  size  h  is  chosen  fine  enough  to  accurately  describe  the 
frequencies  of  interest,  it  may  still  be  many  mean  free  paths  wide.  For  example,  in 
applications  in  radiative  transport,  a  may  be  on  the  order  of  10°  and  the  number  of 
mesh  points  may  be  103.  Then  h  is  on  the  order  of  103  mean  free  paths  wide.  This 
yields  a  discrete  equation  in  which  well  over  half  of  the  singular  values  of  Kh  are 
very  near  unity  and,  because  of  the  discretization  error,  the  smallest  singular  values 
of  Kh  are  well  away  from  unity.  Thus,  even  if  the  lowest  half  of  the  frequencies 
are  resolved  on  a  coarser  grid,  the  effective  condition  number  of  I  —  Kh  due  to  the 
remaining  singular  values  is  extremely  large. 

Simple  source  iteration  is  not  adequate  for  eliminating  the  high-frequency  er¬ 
rors  that  cannot  be  resolved  on  a  coarse  mesh.  Certain  block  relaxation  techniques 
have  recently  been  developed  BOl  that  effectively  attenuate  high-frequency  errors 
regardless  of  the  size  of  h.  Furthermore,  they  have  demonstrated  that  these  re¬ 
laxation  schemes  are  effective  regardless  of  the  anisotropy  of  the  scattering.  Thus 
multigrid  methods  are  more  goner. thy  applicable  than  the  DSA  method. 

Recently,  we  have  developed  .t  relaxation  scheme  on  the  (x,/z)  grid  that  corre¬ 
sponds  to  line  relaxation  in  p  a  triable  shift  [MA6],  We  have  shown  that  this 
yields  a  multigrid  smoothing  rate  ■>:'  1  3  on  the  relevant  singular  vectors.  Moreover, 
after  relaxation  the  error  is  neari>  mbependent  of  /z.  Thus,  a  coarse  grid  correction 
need  not  involve  /z.  More  precisely,  an  effective  coarse  grid  error  correction  equa¬ 
tion  can  be  found  using  the  obser-.  on  that  errors  that  remain  after  relaxation 
are  nearly  independent  of  /z.  A  -<  --.e  based  on  these  observations  is  currently 

being  implemented  to  prove  its  Hf>  -  ss.  The  scheme  has  been  adapted  to  both 

upwind  difference  approximations  ami  the  more  complicated  linear  discontinuous 
difference  schemes. 

Future  Research.  Most  importantly,  the  scheme  can  be  generalized  to  anisotropic 
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scattering  and  to  higher  dimensions.  For  isotropic  scattering,  shifted  line  relaxation 
in  n  requires  solving  diagnonal  block  matrices  perturbed  by  a  rank  one  operator. 
This  can  be  solved  effectively  by  capacitance  matrix  methods  (cf.  [BUlj).  In  the 
anisotropic  case,  the  same  relaxation  scheme  involves  solving  a  simple  block  system 
perturbed  by  a  small  number  of  rank  one  operators.  Although  slightly  more  compli¬ 
cated,  these  may  also  be  solved  by  capacitance  matrix  methods.  In  two  dimensions, 
line  relaxation  in  fj,  becomes  plane  relaxation  in  two  angle  parameters 
Again  these  involve  only  simple  block  matrices  plus  low  rank  perturbations. 

The  algorithm  is  simple  and  inherently  parallel.  First,  the  relaxation  scheme  is 
local  in  space.  Also,  if  a  Jacobi-like  scheme  is  used,  all  lines  can  be  relaxed  simul¬ 
taneously.  If  a  Gauss-Siedel-type  scheme  is  used  where  parallelism  can  be  achieved 
by  using  a  zebra-like  ordering,  every  other  line  can  be  relaxed  simultaneously.  In 
addition,  the  coarse  grids  are  essentially  conventional  diffusion  equations  which  can 
be  handled  by  standard  parallel  multigrid  algorithms. 

The  ultimate  goal  of  this  project  is  the  production  of  software  for  the  solution 
of  these  problems.  The  successful  algorithms  will  be  implemented  in  two  forms. 
One  will  be  designed  for  the  Cray  X  MP  and  directed  toward  production  problems 
at  LANL.  The  second  will  be  coded  for  the  Connection  Machine  (CM2). 

4.  SUPRACONVERGENCE. 

The  classical  approach  to  proving  order  of  accuracy  of  difference  schemes  for 
differential  equations  is  to  examine  the  truncation  error.  If  the  difference  scheme  is 
stable  then  the  order  of  accuracy  is  bounded  by  the  order  of  the  truncation  error. 
Unfortunately,  on  highly  irregular  meshes  the  truncation  error  associated  with  many 
difference  schemes  is  of  lower  order.  For  example,  the  second-divided  difference  has 
second-order  truncation  error  on  a  uniform  mesh,  but  only  first-order  truncation 
error  on  highly  irregular  meshes. 

Much  work  has  been  done  to  alleviate  this  apparent  disadvantage  of  irregular 
meshes.  One  approach  is  to  assume  mat  the  mesh  satisfies  some  smooth  mesh 
function.  Another  approach  is  to  use  more  complicated  difference  schemes,  for 
example,  operator  implicit  schemes. 

Recent  work  by  Manteuffel  and  White  :MAl],  [MA2]  and  Kreiss,  Manteuffel, 
Swartz,  Wendroff  and  White  'KRl  shows  that  for  many  standard  difference  schemes 
the  error  is  actually  second-order  accurate  despite  first-order  truncation  error.  More 
precisely,  consider  the  simple  equa^m 

•/ 


together  with  appropriate  bourn:. .  •  •  r.ditions.  Let  D  be  the  second  divided  dif¬ 
ference  operator  on  the  grid  fur.<  •  •  d  consider  the  discrete  equation 

1>'J  /- 

The  associated  error  equations  can  be  written 

De  -  t 
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where  e  is  the  error  and  t  is  the  truncation  error.  For  irregular  meshes,  jUljoo  =  0(<5) 
where  8  is  the  maximum  mesh  size.  We  know  that 

i;i;u  =  o-'t  x  <  d-'  M  t  ,, . 

Since  \D~  1  ^  is  bounded  independently  of  8,  we  may  conclude  that  e:'x  =  0(8) 
as  well.  However,  the  inequality  above  is  unnecessarily  crude.  It  can  actually  be 
shown  that  ' e\\x  =  0(<52)  despite  the  fact  that  ^  =  0(d).  This  phenomena  has 
been  labeled  supraconvergence  KRl  . 

It  is  clear  that  the  standard  stability/consistency  proof  is  insufficient  to  detect 
supraconvergence.  New  tools  have  been  developed  that  have  led  to  the  discovery 
that  many  schemes  for  ordinary  differential  equations  and  partial  differential  equa¬ 
tions  possess  this  property  on  highly  irregular  meshes.  One  important  tool  is  to 
write  the  differential  equation  as  a  first-order  system  of  equations  using  auxiliary 
unknowns,  discretize  the  first-order  system  with  simple  difference  schemes  and  then 
algebraically  eliminate  the  auxiliary  variables.  This  process,  if  done  correctly,  yields 
a  compact-as-possible  difference  scheme  on  the  original  equation.  Since  the  solution 
to  the  reduced  equations  is  algebraically  identical  to  the  solution  of  the  first-order 
system,  it  inherits  any  properties  that  the  first-order  system  possesses. 

Results.  The  theory  of  supraconvergence  has  been  applied  to  higher-order  ordi¬ 
nary  differential  equations  with  surprising  results.  It  is  quite  simple  to  write  a 
higher-order  equation  as  a  first  order  system  and  to  then  discretize  this  system  us¬ 
ing  a  second  order  accurate  difference  scheme.  In  [MA3j  (attached),  a  calculus  of 
difference  operators  is  introduced.  Auxiliary  operators  are  constructed  that  allow 
for  easy  reduction  of  the  discrete  system  to  a  discrete  compact-as-possible  system 
for  the  original  equation.  A  more  detailed  paper  [MA4j  which  is  near  completion 
shows  that  this  process  is  essentially  unique.  The  process  yields  second-order  ac¬ 
curate  difference  schemes  for  any  order  ODE.  The  truncation  error,  however,  is 
not  second  order  on  irregular  meshes.  In  fact,  for  say  a  fourth  order  equation,  the 
truncation  error  is  0(8~2),  not  only  inconsistant  but  increasing  as  8  becomes  small. 

Boundary  conditions  for  partial  differential  equations  become  especially  trans¬ 
parent  in  the  context  of  first-order  systems.  The  algebraic  reduction  of  the  first- 
order  system  yields  correct  boundary  conditions  for  the  original  equation.  This  im¬ 
portant  aspect  of  the  research  will  have  implications  for  choosing  difference  schemes 
to  link  composite  grids. 

Future  Research.  These  same  tools  have  produced  results  for  partial  differential 
equations.  On  tensor  product  meshes  the  extension  is  easy  [MA2].  However,  on 
more  complicated  meshes  the  problem  is  difficult. 

A  recent  report  (LElj  develops  a  r  lass  of  compact-as-possible  difference  schemes 
for  parabolic  equations  on  time  lines,  but  irregular  spacial  mesh.  The  class  possesses 
second-order  accuracy  despite  first-order  truncation  errors.  This  report  also  includes 
results  on  hyperbolic  equations.  We  -how  that  a  class  of  conservative  difference 
schemes  are  first-order  accurate  on  irregular  meshes  despite  inconsistent  truncation 
error.  These  results  are  extended  to  the  nonlinear  wave  equation. 
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The  more  compelling  question  of  supraconvergence  for  elliptic  equations  on  ir¬ 
regular  meshes  remains  open.  Current  research  that  also  involves  Andrew  B.  White, 
Jr.,  in  Los  Alamos  National  Laboratory,  Max  Gunzberger  at  Carnegie-Mellon  Uni¬ 
versity  and  David  Levermore  of  the  University  of  Arizona  is  focusing  on  9-point 
schemes  on  logically  rectangular  meshes  and  variable  point  schemes  on  triangular 
meshes.  The  main  idea  is  to  define  auxiliary  unknowns  that  are  not  necessarily 
associated  with  the  nodes  of  the  mesh,  to  form  two  different  inner  product  spaces, 
one  that  includes  the  original  unknowns  and  another  that  includes  the  auxiliary 
unknowns,  and  to  determine  the  adjoints  of  the  difference  schemes  in  the  respective 
inner  products.  If  the  truncation  error  can  be  conveniently  decomposed  into  por¬ 
tions  in  the  range  and  the  null  space  of  the  adjoint,  then  the  truncation  error  can  be 
written  as  a  product  of  the  first-order  system  times  terms  of  possibly  higher-order 
plus  terms  in  the  null  space.  In  this  way  schemes  can  be  analyzed  to  see  if  they 
possess  supraconvergence. 
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